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Details  of  the  numerical  solution  techniques,  including  explicit,  point-implicit,  and  fully 
implicit  schemes,  used  in  a  new  quasi  two-dimensional  procedure  for  the  transient  solution 
of  real  fluids  in  system  lines  and  volumes  are  presented.  The  procedure  is  coupled  with  a 
real-fluid  properties  database  so  that  both  compressible  and  incompressible  fluids  may  be 
considered  using  the  same  code.  The  procedure  has  been  implemented  in  Matlab/Simulink® 
as  well  as  Fortran95  to  allow  for  application  on  a  wide  variety  of  computer  platforms.  The 
computational  efficiency  of  the  various  numerical  methods  is  discussed  to  aid  in  selection  for 
specific  appfications.  Resufts  for  the  transient  flows  of  gaseous  nitrogen  and  water  in  a  pipe 
are  presented  to  demonstrate  the  capability  of  the  current  techniques  and  the  unsteady  flow 
physics  that  can  occur  in  system  lines. 
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Nomenclature 

=  cross-sectional  area 
=  wetted  area 

=  relative  acceleration  between  inertial  and  absolute  frames  of  reference 
=  Courant-Friedrichs-Lewi  number 
=  speed  of  sound 
=  total  energy  =  p  (e  +  if  IF) 

=  internal  energy 

=  flux  of  mass,  momentum,  and  energy 
=  static  enthalpy 
=  cell  index 

=  temporal  order  of  accuracy  for  ODE15s  solution  methodology 
=  index  variables 
=  static  pressure 
=  heat  flux  input 

=  source  vector 
=  time 
=  state  vector 
=  absolute  velocity 
=  relative  velocity 
=  fluid  volume  =  Up 
=  distance  along  solution  domain 
=  fluid  density 
=  shear  force  (friction) 


Professor,  Mechanical  and  Aeronautical  Engineering,  AIAA  Associate  Fellow. 

1  Engineering  Specialist,  Systems  Engineering,  P.O.  Box  13222,  Dept.  5271,  AIAA  Member. 


1 

American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


I.  Introduction 

THE  simulation  of  unsteady  flows  of  non-idealized  fluids  in  system  lines  and  volumes  is  of  interest  for  a  wide 
variety  of  applications  and  industries,  including  air-conditioning  systems,  water/steam/oil  piping  networks, 
refinery  systems,  gas-turbine  secondary  flow-path  and  cooling  networks,  and  liquid  rocket  engine  propellant  lines. 
Several  approaches  for  simulating  the  dynamic  behavior  of  such  fluid-transmission  lines  have  been  reported.  The 
lumped-analysis  approach  treats  a  flow  passage  as  a  series  of  fluid  control  volumes  that  conserve  mass  and  energy 
linked  by  flow  resistance  elements  that  compute  the  flow  between  the  volumes 1 .  While  this  approach  does  conserve 
momentum  in  a  quasi-steady  sense  at  the  flow  resistances,  the  unsteady  momentum  term  in  the  governing  equations 
is  neglected.  Although  so-called  “continuity”  waves  can  be  captured  using  this  approach,  neglecting  these  terms 
leads  to  an  inability  to  capture  true  “dynamic”  waves  required  for  simulating  such  phenomenon  as  water-hammer 
and  pressure  surge2.  Another  approach  utilizes  the  Method  of  Characteristics3,  which  is  a  general  method  for 
solving  particular  types  of  differential  equations.  The  governing  equations  for  fluid  flow  are  compatible  with  this 
method  and  it  has  been  used  for  simulation  of  fluid  transmission  lines4.  While  the  unsteady  momentum  terms  are 
retained  using  this  method,  other  problems,  particularly  at  the  boundaries  of  components,  make  it  difficult  to  apply 
to  a  modular  system-level  simulation  tool.  Modal  methods  have  also  been  used  when  solution  in  the  frequency 
domain  is  possible5.  This  technique  represents  the  pressure  and  velocity  distributions  in  the  flow  domain  as  a  sum 
of  an  infinite  series  of  mode  shapes,  similar  to  a  Fourier  series  solution.  While  this  method  does  present  an  elegant 
and  efficient  method  for  simulating  idealized  flows  (e.g.  incompressible,  inviscid,  laminar,  etc),  the  addition  of 
turbulent  flow,  real-fluid  properties,  heat  transfer  and  phase  change  complicate  the  application  of  the  method  and 
reduce  its  attractiveness. 

Since  unsteady  phenomena  such  as  wave  dynamics  play  an  important  role  in  the  operation  and  testing  of 
systems  that  contain  fluid  lines,  a  method  that  captures  these  transient  effects  is  required.  The  development  of  a 
quasi  two-dimensional,  unsteady,  two-phase  flow  solver  with  heat  transfer  and  real-fluid  properties  using  standard 
finite-difference/control-volume  solution  methods  has  been  the  subject  of  the  present  effort.  Such  a  solver  is  geared 
towards  modeling  the  dynamic  behavior  of  fluid-filled  lines  and  passages  (i.e.  the  solution  domain  is  much  larger  in 
one  spatial  dimension  than  in  the  others)  accounting  for  the  effects  of  changing  cross-sectional  area.  In  addition,  the 
solver  must  be  suitable  for  use  as  a  module  in  larger  system-level  transient  simulations  of  hydraulic  and  pneumatic 
systems,  so  the  solution  method  must  be  computationally  efficient.  The  following  sections  describe  the  modeling 
approach,  numerical  methodologies  and  test  cases  that  have  been  utilized  during  the  development  of  this  model. 
Results  are  then  shown  for  transient  pipe  flow  of  both  nitrogen  and  water  as  a  demonstration  of  the  numerical 
capability  and  fluid  physics  that  can  be  captured  with  the  current  procedure. 

II.  Approach 

The  model  developed  here  represents  fluid  lines  and  flow  passages  where  the  length  of  the  domain  is  much 
larger  than  the  hydraulic  diameter  of  the  domain  as  a  quasi  two-dimensional  domain.  For  these  types  of 
components,  flow  separations  and  non-axial  velocities  are  minimal,  hence  the  quasi-2D  assumption  is  valid.  The 
solver  is  targeted  to  the  commercial  Simulink"  dynamic  simulation  software  package  from  the  MathWorks  for 
integration  into  a  larger  suite  of  modules  developed  for  simulating  various  systems.  Simulink1'  was  selected  since  it 
offers  a  wide  range  of  capabilities,  over  a  dozen  robust  differential  equation  solvers,  extensive  documentation  and 
technical  support,  a  modern  graphically-based  modeling  paradigm,  an  existing  user  community  across  many 
disciplines,  and  commercially- funded  code  development  and  maintenance.  A  Fortran95  code  using  more  traditional 
solution  methods  is  also  being  developed  in  parallel  to  provide  test  cases  for  the  Simulink"  module. 

The  solver  is  being  developed  to  account  for  varying  flow  area,  friction,  losses  due  to  bends,  real-fluid  and  two- 
phase  flow  effects,  gravity  and  acceleration,  heat  transfer,  and  the  capability  to  produce  unsteady  and  steady-state 
solutions.  Fluid  properties  are  obtained  from  the  REFPROP  fluid  property  database6  available  from  the  National 
Institute  for  Standards  and  Testing  (NIST).  This  database  utilizes  state-of-the-art  Equation  of  State  models  to  fully 
describe  “real  fluid”  properties  over  a  wide  range  of  thermodynamic  conditions,  including  liquid,  vapor,  mixed 
phase  and  supercritical  fluid  regimes.  Properties  that  completely  define  the  fluid  thermodynamic  state,  as  well  as 
transport  properties,  are  available  as  a  function  of  any  two  thermodynamic  parameters.  Validated  fluid  models  for 
over  80  pure  fluids  and  over  180  fluid  mixtures  are  available  in  the  database.  The  database  is  accessed  through  a 
suite  of  Fortran77  subroutines  that  were  compiled  into  the  codes  and  obtain  fluid  equation  of  state  model  parameters 
from  external  files  provided  with  the  database. 

Friction  and  heat  transfer  are  modeled  as  source  terms  in  the  governing  equations.  This  approach  allows  the 
flow  to  be  modeled  as  one-dimensional  and  facilitates  the  computational  efficiency.  Friction  (i.e.  viscous  losses  and 
losses  due  to  bends  and  fittings)  and  heat  transfer  coefficients  are  obtained  from  suitable  correlations  between  the 
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flow  variables  and  the  source  terms.  The  methodology  for  modeling  two-phase  flow  will  be  included  in  the  final 
paper. 


III.  Governing  Equations 

The  governing  equations  consist  of  the  quasi  two-dimensional  Euler  equations  with  source  terms  for  viscous 
effects  and  general  acceleration  in  a  non-inertial  frame. 
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where  p  is  the  density,  u  is  the  absolute  velocity  in  the  inertial  frame,  u,  is  the  relative  velocity  in  the  non-inertial 
frame,  ar  is  the  gravitational  and  relative  acceleration  in  the  non-inertial  frame,  E  is  the  total  energy,  p  is  the 
pressure,  H  is  the  stagnation  enthalpy,  V  is  the  local  cell  volume,  A  is  the  local  line  cross-sectional  area,  dx  is  the  cell 
line  length,  Ap  is  the  perimeter  area,  xu  is  the  wall  shear  stress.  The  wall  shear  stress  can  be  written  in  terms  of  a 
friction  factor,/  which  is  a  function  of  the  local  Reynolds  number  and  the  wall  surface  roughness.  For  the  viscous 
flow  examples  presented  below,  the  Churchill  correlation7  was  used  to  determine  single-phase  friction  factors  using 
the  equation: 
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where  Re  is  the  local  Reynolds  number  based  on  diameter. 

Closure  between  the  pressure,  density,  and  energy  is  performed  with  an  equation  of  state  for  real  fluids  by  using 
the  NIST  REFPROP  suite  of  thermodynamic  routines.  This  allows  for  the  solution  of  compressible,  incompressible, 
and  two-phase  flows. 


IV.  Numerical  Methods 

Equations  (1)  and  (2)  may  be  solved  with  a  variety  of  numerical  schemes  that  are  second-order  accurate  in 
space  and  time  for  “steady”  or  transient  flow.  These  schemes  currently  include  implicit,  point-implicit,  and  explicit 
time-marching  procedures.  Multiple-grid  convergence  acceleration  may  be  used  in  steady-flow  schemes  as  well  as 
in  the  inner  iteration  of  time-accurate  point-implicit  schemes  to  dramatically  reduce  the  number  of  steady  or  inner 
iterations.  Multiple  numerical  schemes  have  been  developed  in  the  current  procedure  to  provide  flexibility  and 
robustness  when  solving  different  flow  conditions. 

A.  Implicit 

The  fully  implicit  discretization  of  Eqs.  (1)  and  (2)  may  be  achieved  by  representing  the  primary  variables 
contained  in  U  at  the  new  time  step  rather  than  at  the  current  time  step.  This  results  in  an  N  system  of  equations 
where  N  is  the  number  of  nodes  in  the  line  to  be  solved.  The  coefficients  of  the  matrix  which  multiplies  the  primary 
variable  column  matrix  can  be  written  as  a  tri-diagonal  system  (for  second  order  spatial  accuracy).  Solution  of  the 
system  of  equations  can  then  be  determined  by  finding  the  inverse  of  the  coefficients  matrix. 

An  alternative  to  solving  the  system  fully  coupled  in  terms  of  a  block-tridiagonal  system,  Eqs.  (1)  and  (2)  may  be 
solved  iteratively  in  a  sequential  manner.  The  Matlab/Simulink  solution  procedure8,  ODE15s,  treats  the  discretized 
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form  of  Eqs.  (1)  and  (2)  as  a  sequential  set  of  ordinary  differential  equations  and  solves  them  iteratively  with  a 
Newton  iteration  according  to  the  update  formula: 


-At 


* +s 

dx 


=  0 
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where  k  is  the  temporal  order  of  accuracy.  In  the  current  investigation,  both  the  temporal  and  spatial  accuracy  is 
second  order.  The  time-step  size  is  determined  by  the  Courant-Friedrichs-Lewi  condition  i.e. 
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where  CFL  is  the  stability  number  and  c  is  the  local  speed  of  sound.  The  CFL  stability  number  for  the  implicit 
technique  can  be  very  large.  Typical  values  are  100-10,000. 

B.  Explicit 

The  explicit  time-marching  of  Eqs.  (1)  and  (2)  may  also  be  performed  using  a  Lax-Wendroff  scheme,  as 
described  by  Ni9.  In  this  procedure,  time-marching  the  primary  variables,  U,  corresponds  to  an  update  in  time  at 
each  node  according  to  a  second  order  temporally  accurate  Taylor  series  formula: 
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The  Ni  scheme  is  a  finite  volume  integration  method.  The  scheme  is  applied  in  a  three  step  process.  First  the  change 
is  at  the  center  of  each  computational  cell  is  approximated  as 
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The  correction  to  the  points  i  and  i+1  are  determined  using  the  distribution  formulas  given  below: 


w>4 


AU-—AF 

Ax 


AU  +—AF 
Ax 


where 


A F  = 


f  dF^ 


\dU  j 


AU 


The  primary  variables  are  then  updated  by 

dU^idUXu+idUiJ^, 


8Un+l=Un+dU 


(8) 

(9) 


(10) 


(11) 

(12) 


4 

American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


The  first  order  changes  in  time,  dU,  are  found  directly  from  application  of  Eqs.  (1)  and  (2)  on  the  primary 
computational  cells.  The  second  order  changes  in  time,  oU,  represented  by  the  second-order  derivatives  of  U 
above,  are  found  by  applying  the  time-derivative  of  Eqs.  (1)  and  (2)  on  computational  cells  that  are  translated  'A  cell 
line  length  from  the  primary  computational  cells.  The  CFL  stability  number  that  sets  the  time-step  size  is  set  to  0.7 
due  to  the  explicit  nature  of  this  scheme. 

For  steady  flows  or  the  inner  iteration  of  the  point-implicit  technique  described  below,  a  multiple-grid 
acceleration  procedure  described  by  Ni  is  used  to  greatly  reduce  the  number  of  time-steps  to  reach  convergence. 
This  allows  the  scheme  to  be  competitive  with  the  implicit  procedure  described  above. 

C.  Point  Implicit 

A  point-implicit  procedure  can  be  made  up  from  either  the  implicit  or  explicit  schemes  described  above.  For  this 
procedure,  the  time-rate  changes  in  Eqs.  (4)  or  (6)  are  assumed  to  be  in  pseudo-time.  A  true  time  derivative  of  the 
primary  variables  is  added  to  the  right-hand  side  of  Eq.  (1).  This  true  time  derivative  is  derived  from  stored 
solutions  of  the  primary  variables  at  k+1  previous  time-steps,  where  again,  k  is  the  temporal  order  of  accuracy 
desired.  An  inner-iteration  is  then  used  to  drive  the  right-hand-side  to  zero.  The  advantage  of  this  approach  is  that 
steady-flow  acceleration  techniques  may  be  used  during  the  inner  iteration  such  as  a  multiple-grid  scheme.  Also, 
this  technique  provides  flexibility  in  performing  multiple  lines  of  a  network  in  parallel.  This  procedure  is  known  as 
the  dual-time-step  technique10  and  has  many  similarities  to  the  Newton  iteration  used  in  the  implicit  procedure. 

V.  Results 

Transient  pipe  flow  simulations  of  gaseous  nitrogen  and  liquid  water  have  been  performed  to  demonstrate  the 
capability  of  the  present  procedure.  A  0.254  m  (10  in)  pipe  with  a  0.0254  m  (1  in)  diameter  was  used  for  these 
simulations.  The  surface  roughness  of  the  pipe,  necessary  for  viscous  flow  simulations,  was  2.54  x  10"6  m 
(0.0001  in).  The  inlet  total  pressure  was  held  at  3.45  kPa  (50  psia)  as  well  as  the  initial  exit  static  pressure.  The 
inlet  total  temperature  was  held  at  300°K  (540°R).  Thus,  the  initial  velocity  of  the  flow  in  the  pipe  was  zero.  At 
time  0.0001  seconds,  the  exit  pressure  was  reduced  to  3.31  kPa  (48  psia).  This  sudden  reduction  in  pressure  initiates 
a  series  of  expansion  and  compression  waves  through  the  pipe  and  an  increase  in  the  velocity  of  the  flow.  For 
inviscid-flow  cases,  the  pressure  decays  in  the  pipe  to  a  uniform  value  and  the  velocity  grows  asymptotically  to  a 
uniform  value.  For  viscous-flow  cases,  the  pressure  in  the  pipe  decays  but  to  a  non-uniform  distribution  with  the 
inlet  pressure  higher  than  the  exit  due  to  friction  losses.  The  velocity  of  the  pipe  flow,  however,  will  once  again 
grow  asymptotically  until  it  reaches  a  uniform  value.  All  simulations  were  run  in  time  until  both  the  pressure  and 
velocity  distributions  through  the  pipe  remained  time-independent.  A  total  of  65  grid  points  with  uniform  spacing 
was  used  along  the  pipe  in  order  to  discretize  the  governing  equations. 

A.  Transient  Nitrogen  Pipe  Flow 

The  pipe  flow  simulation  with  nitrogen  is  meant  to  test  the  current  procedure’s  capabilities  for  compressible 
flow.  The  explicit  and  point-implicit  numerical  procedures  were  both  run  in  order  to  make  comparisons  and  to 
demonstrate  features  of  each  technique.  Figure  1  shows  the  predicted  pressure  and  velocity  as  a  function  of  time 
using  the  explicit  numerical  technique  for  inviscid  flow  at  five  locations  along  the  pipe.  In  this  scheme,  the  time 
step  size  is  determined  by  the  minimum  time  step  of  all  the  cells  in  the  domain.  The  sinusoidal  variation  of  pressure 
with  time  is  evidence  of  the  expansion  and  compression  waves  propagating  through  the  pipe.  The  time-variation  of 
the  velocity  corresponds  to  these  waves.  The  pressure  decays  in  time  until  it  reaches  the  exit  pressure  throughout 
the  pipe.  The  velocity  increases  until  it  reaches  a  value  of  approximately  31  m/s  (102  ft/s). 


5 

American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


a)  Pressure 

Figure  1  -  Nitrogen  Transient  Inviscid  Pipe  Flow  using  Explicit  Numerical  Technique 


Figure  2  shows  the  predicted  results  using  the  point-implicit  numerical  procedure.  In  this  procedure,  the  time- 
step  size  can  be  prescribed  by  the  user  to  be  any  value  that  results  in  a  desired  solution  time  and  resolution  of 
temporal  features  in  the  flow.  For  the  present  case,  the  time-step  size  was  set  at  0.01  seconds  which  was 
approximately  1,500  times  larger  than  that  required  for  in  the  explicit  procedure.  The  capability  to  use  a  large  time- 
step  size  in  the  point-implicit  procedure  allows  for  very  fast  solution  time.  Flowever,  as  shown  in  Figure  2, 
resolution  of  the  high-frequency  waves  and  their  effect  on  the  pressure  and  velocity  is  lost  as  a  result  of  the 
increased  time-step  size. 


Nitrogen,  Inviscid,  Pt  Implicit 
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a)  Pressure 

Figure  2  -  Nitrogen  Transient  Inviscid  Pipe 


Nitrogen,  Inviscid,  Pt  Implicit 
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b)  Velocity 

using  Point-Implicit  Numerical  Technique 


Both  numerical  schemes  were  run  with  viscous  shear  stress  effects  modeled  with  a  friction  factor  correlation  as 
described  above  as  a  further  demonstration  of  the  current  procedure’s  capability.  Figure  3  shows  the  predicted 
pressure  and  velocity  resulting  from  the  explicit  transient  viscous  flow  simulation  of  nitrogen.  The  pressure  decays 
to  a  non-uniform  distribution  through  the  pipe  for  viscous  flow.  The  pressure  drop  along  the  pipe  agrees  very  well 
with  what  is  expected  from  a  lumped  parameter  analysis  with  the  energy  equation.  The  velocity  asymptotically 
increases  with  time  to  a  level  below  that  predicted  without  viscous  effects  (shown  in  Figure  1). 
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a)  Pressure 

Figure  3  -  Nitrogen  Transient  Viscous  Pipe  Flow  using  Explicit  Numerical  Technique 


The  nitrogen  viscous  flow  results  using  the  point-implicit  numerical  procedure  are  shown  in  Figure  4.  The 
transient  solution  is  obtained  with  a  time-step  size  1,500  times  greater  than  that  for  the  explicit  solution  shown  in 
Figure  3  but  with  the  loss  in  temporal  resolution  of  the  pressure  and  velocity  solution. 


a)  Pressure  b)  Velocity 

Figure  4  -  Nitrogen  Transient  Viscous  Pipe  Flow  using  Point-Implicit  Numerical  Technique 

Transient  Water  Pipe  Flow 

A  similar  set  of  simulations  were  performed  using  the  present  procedure  with  water  to  demonstrate  the  capability 
to  predict  transient  flows  in  incompressible  fluids.  The  main  critical  difference  with  incompressible  fluids  is  that  the 
speed  of  sound  and  wave  speeds  are  much  larger  for  incompressible  fluids.  This  makes  the  use  of  explicit  numerical 
techniques  time-consuming.  Resolution  of  the  unsteady  sinusoidal  transients  shown  above  are  possible,  but  quite 
time  consuming.  The  effect  of  the  high-frequency  waves  are  not  always  of  interest  so  that  numerical  techniques  that 
can  take  large  time-steps  and  reduce  solution  time  are  often  attractable. 

Figure  5  shows  the  transient  point-implicit  inviscid  solution  of  water  through  the  previously  described 
configuration.  Note  that  the  time  to  reach  equilibrium  for  water  is  much  longer  than  that  required  for  nitrogen  due  to 
the  incompressibility  of  the  fluid. 
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a)  Pressure  b)  Velocity 

Figure  5  -  Water  Transient  Inviscid  Pipe  Flow  using  Point-Implicit  Numerical  Technique 

The  size  of  the  time  step  used  in  the  point-implicit  solution  was  0.01  seconds.  The  time  step  size  required  in  the 
corresponding  explicit  solution  procedure  was  approximately  2xl0’6  seconds  or  5,000  times  smaller  than  that  used  in 
the  point-implicit  solution.  As  a  result,  over  600,000  time  steps  would  be  required  to  solve  the  1 .2  seconds  of  time 
to  reach  a  time-independent  solution. 

The  transient  viscous  solution  for  water  pipe  flow  using  the  point-implicit  procedure  is  shown  in  Figure  6.  The 
viscous  pressure  solution  converges  to  a  non-uniform  time-independent  distribution  that  agrees  very  well  with 
lumped  parameter  energy  equation  analysis.  The  velocity  converges  to  a  uniform  time-independent  solution  with  a 
slightly  lower  velocity  than  the  predicted  inviscid  time-independent  solution. 


Figure  6  -  Water  Transient  Viscous  Pipe  Flow  using  Point-Implicit  Numerical  Technique 

As  for  the  inviscid  case,  the  time  step  size  required  for  the  explicit  solution  is  approximately  5,000  times  smaller 
requiring  an  excessive  solution  time.  Both  inviscid  and  viscous  solutions  have  been  run,  however,  to  verify  solution 
capability.  Comparison  between  the  explicit  and  point-implicit  solutions  are  similar  to  that  for  nitrogen  with  the 
explicit  solutions  showing  very  high  frequency  content  due  to  the  various  expansion  and  compression  waves 
propagating  through  the  pipe. 


VI.  Conclusion 

A  new  baseline  procedure  has  been  developed  for  the  numerical  solution  of  transient  quasi  two-dimensional  flow 
in  engine  and  system  lines  and  volumes.  This  new  procedure  is  driven  written  in  both  Matlab/Simulink  and 
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Fortran95.  A  variety  of  numerical  solution  algorithms  are  included  to  allow  for  various  needs  depending  on  the 
desired  solution  time  and  frequency  resolution.  Currently,  explicit  and  point-implicit  numerical  techniques  have 
been  implemented  and  validated.  The  point-implict  technique  is  very  useful  for  finding  rapid  turn-around  of 
solutions  at  the  expense  of  some  loss  in  temporal  resolution.  The  explicit  technique  is  very  useful  for  predicting  the 
high-frequency  content  of  transient  flow  fields.  Implementation  of  a  fully  implicit  numerical  technique  is  underway 
and  results  will  also  be  shown  for  this  technique  in  the  final  paper. 
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